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ABSTRACT 

In this paper we present a characteristic method for solving the transfer equation 
in differentially moving media in a curved spacetime. The method is completely 
general, but its capabilities are exploited at best in presence of symmetries, when 
the existence of conserved quantities allows to derive analytical expressions for 
the photon trajectories in phase space. In spherically-symmetric, stationary 
configurations the solution of the transfer problem is reduced to the integration 
of a single ordinary differential equation along the bi-parametric family of 
characteristic rays. Accurate expressions for the radiative processes relevant to 
continuum transfer in a hot astrophysical plasma have been used in evaluating 
the source term, including relativistic e-p, e-e bremsstrahlung and Compton 
scattering. A numerical code for the solution of the transfer problem in moving 
media in a Schwarzschild spacetime has been developed and tested. Some 
applications, concerning "hot" and "cold" accretion onto non-rotating black 
holes as well as static atmospheres around neutron stars, are presented and 
discussed. 

Subject headings: accretion, accretion disks - numerical methods - radiative 
transfer - relativity 
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1. INTRODUCTION 

Radiative transfer in high energy, fast moving plasmas in a strong gravitational 
field is today at the basis of a large number of currently interesting astrophysical 
applications; accretion onto compact objects, jets, stellar collapse and supernova 
expanding envelopes are just some examples of this. 

Since the pioneering works by Thomas (1930), Simon (1963) and Lindquist 
(1966), astrophysical relativistic transfer received wide attention (see e.g. 
Mihalas & Mihalas 1984 for references to earlier papers). It was realized 
long ago (see Castor 1972, Mihalas 1980 and references therein) that, for 
relativistic flows, the interaction between matter and radiation is most easily 
described if the material properties and the radiation field are evaluated in 
the frame in which the medium is at rest. The comoving frame transfer 
equation (CTE) has been considered by Mihalas (1980), Hauschildt & Wehrse 
(1991) in the framework of special relativity, and by Schmidt-Burgk (1978), 
Thorne (1981), Schinder & Bludman (1989) in the general-relativistic case. 
Different approaches for the solution of the relativistic transfer problem in 
planar or spherical geometry have been suggested. They can be grouped, 
schematically, into three wide classes: direct solution of the CTE using 
discretization techniques, moment expansion and integration of the CTE along 
characteristic directions. 

The solution of the CTE by finite differencing, like in the DOME method 
(Hauschildt & Wehrse 1991), works well in geometrically thin layers, but 
the treatment of extended atmospheres requires a prohibitive number of 
discrete elements to obtain a fair angular resolution. In the relatively simple 
case examined by Hauschildt & Wehrse, the numerical calculations must be 
performed on supercomputers even for low resolution grids. 

The expansion of the specific intensity in spherical harmonics (moments) 
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has the main advantage of reducing the dimensionahty of the problem since 
the angular dependence is suppressed. On the other hand the solution of the 
transfer problem is reconduced to the solution of a recursive system of partial 
differential equations that must be truncated at a given order, introducing a 
certain number of closure conditions. This approach is at the basis of the 
flux limited diffusion theory (FDT) developed by Levermore & Pomraining 
(1981) and generalized by Pomraning (1983), Anile & Sammartino (1989) 
and Anile & Romano (1992). Although in the gray case the FDT provides 
a self-consistent closure function by solving the differential equation for the 
flux limiter, the extension to the frequency-dependent problem seems far from 
being obvious. A very sophisticated, general-relativistic version of the moment 
formalism was presented by Thorne (1981). It is based on an expansion in 
projected, symmetric, trace-free (PSTF) moments and, upon truncation, the 
resulting system of equations can be solved introducing the required number of 
closure conditions. The closure functions must be specifled "a priori" and should 
reproduce the correct asymptotic limits for the radiation moments when free 
streaming and diffusion are approached. This method has been fruitfully applied 
to the solution of astrophysical problems with planar or spherical symmetry, 
both in the gray and in the frequency-dependent case (see e.g. TuroUa & 
Nobih 1988, Nobih, TuroUa & Zampieri 1993, Zampieri, TuroUa & Treves 1993). 
While the arbitrariness of the closure functions is not a serious problem in 
the gray case, where a large number of moments can be used, it becomes a 
major complication when frequency-dependent transfer is tackled. In fact, to 
make the numerical solution affordable, only the flrst two moments can be 
taken into account, so that the choice of the closures has a non negligible 
impact on the results. Moreover, the extension to the bidimensional case is 
very complicated and requires the speciflcation of 12 closure relations, making 
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the method unacceptably dependent on the choice of a large number of free 
functions. 

Characteristics methods are based on the fact that the transfer equation 
is just the Boltzmann equation for the photon distribution function in phase- 
space (see e.g. Lindquist 1966). The hyperbohc character of the Boltzmann 
equation implies that the CTE can be always reduced to a single ordinary 
differential equation along the characteristic rays. The tangent ray method 
(TRM) developed in a series of papers by Mihalas and coworkers (Mihalas, 
Kunasz & Hummer 1975, 1976a, b, Mihalas 1980) uses a semi-characteristic 
approach in which the integration is performed along the characteristics of 
the "spatial" part of the differential operator (the tangent rays), while the 
frequency derivative is treated by means of a standard finite-differences scheme. 
A fully characteristic method for the solution of the general-relativistic transfer 
problem has been discussed by Schmidt-Burgk (1978), Schinder (1988) and 
Schinder & Bludman (1989). All these investigations dealt with stationary, 
spherically symmetric space-times, which admit three Killing vectors: the 
existence of the associated constants of the motion can be used to obtain 
simple expressions for the characteristic rays. The analysis by Schinder & 
Bludman (1989) was actually restricted to a spacetime characterized by a 
stationary lagrangian line-element, which corresponds to a vanishing eulerian 
velocity field for the matter configuration; their test models refer, in fact, to a 
static atmosphere. The work by Schmidt-Burgk (1978), although finalized to 
accretion onto a Schwarzschild hole, is, to our knowledge, the only example of an 
exact solution of the CTE taking into account both the effects of dynamics and 
strong gravity. For their simple mathematical structure, characteristic methods 
seem to be promising to cope with realistic astrophysical problems. Moreover, 
they can be quite naturally extended to more than one spatial dimension, the 
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major complication coming from the higher number of ODEs that must be 
solved to compute the characteristic trajectories. 

Previous investigations were mainly concerned with the development of 
efficient methods for the solution of the CTE, assuming rather simple, often "ad 
hoc", expressions for the emission and absorption coefficients. This approach 
is completely justified if one is interested in investigating the formation of 
particular spectral features, like lines or absorption edges. On the other 
hand, in all situations in which attention is focussed on the continuum, 
an accurate treatment of all relevant radiative processes becomes important. 
When dealing with hot plasmas, the dominant radiative processes are non- 
conservative scattering and bremsstrahlung. Solutions presented by Schmidt- 
Burgk (1978) refers to a hot, magnetized plasma and takes into account 
scattering absorption and synchrotron absorption/emission; the coUisional term 
in the Boltzmann equation is written using suitable approximations. On the 
other hand, approaches based on moment expansion (see e.g. Pomraining 1973 
and references therein, Thorne 1981, Prasad et al. 1988) do not permit an 
exact description of anisotropic and non-coherent scattering, usually treated 
in the Fokker-Planck approximation. A rigorous treatments of the Compton 
scattering can be found in Kershaw, Prasad & Season (1986), Kershaw (1987), 
Shestakov, Kershaw & Prasad (1988), but their results are never been included 
in transfer codes devoted to astrophysical applications. 

In the following we discuss a fully characteristic approach to the solution 
of the transfer equation in its more general form; results are then specialized to 
stationary, spherically-symmetric or plane-parallel configurations. Particular 
care will be devoted to a detailed treatment of the source term for an 
unmagnetized, fully ionized, non-degenerate hydrogen gas. A numerical code 
is described and applications to accretion onto black holes and neutron stars 
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are finally presented. 



2. RADIATIVE TRANSFER 

In this section we consider the characteristic form of the radiative transfer 
equation in the comoving frame. In subsection a) the CTE and the equations 
for its characteristic trajectories are derived in the more general case, when 
no symmetries are present. A particularization to spherically-symmetric 
configurations is presented in subsection b) . Finally in subsection c) the choice 
of boundary conditions is discussed. In relativistic transfer the radiation field 
is naturally described by the photon distribution function in the phase-space, 
/, that is related to the specific intensity by 2/ = c^I/h^v^. Geometrized 
units {c = G = h = 1) are used throughout and lengths are in units of the 
gravitational radius r = 2M. 



The relativistic transfer equation, written in covariant form, is just the 
Boltzmann equation for /(x, p) 



where p = dx/d^ is the photon 4-momentum, ^ is an affine parameter along the 
null geodesic and the coUisional term g accounts for the interactions between 
matter and radiation (see e.g. Lindquist 1966, Thorne 1981). The differential 
operator in equation (1) acts not merely in spacetime but in the full photons 
phase-space, made up by the spacetime plus the null tangent space at each 
point along the photon trajectory. 

Since /(x, p) is a relativistic invariant, equation (1) holds in any frame. 
However, the material properties (e.g. opacity and emission coefficients, 
scattering cross-section), which enter the expression of the source term are 



a) The Radiative Transfer Equation 



^(x, p) 
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naturally defined with respect to observers who are locally and instantaneously 
at rest with the matter (LRF). In the following we adopt a fiducial observer 
comoving with the fiuid, which carries a tetrad ea and has 4-velocity u = 
eg. If spacetime, matter and radiation share some common symmetries, the 
orientation of some of the spatial vectors of the tetrad follows in a natural way. 
For example, in spherical symmetry, as it will be discussed in detail later on, it 
is convenient to chose ej orthogonal to the 9 and coordinate directions. With 
respect to the tetrad, the components of the photon 4-momentum are 

/ = {E,Eii,E{l- ii'^Y/'^cos^,E{\- li^f/'^si-n^) (2) 

where E is the photon energy, ji is the cosine of the angle between the photon 
direction and ej, and $ is the corresponding azimuthal angle, all measured in 
the LRF. The three quantities jj, and $ have an immediate physical meaning 
and they will be used as independent variables (momentum variables) together 
with the spacetime coordinates a;* to tick events on the light-cone of the phase- 
space. The total derivative in equation (1) can be explicitated as 



df df dp' 
iP + 



9.. 9pa ^ 

df i df dE df dn df d$ _ 
^ dE~di ^ djld^ ^ d^~d^ ~ ^ ' 

where = p"e^. The variation of E, fj. and $ along the photon trajectory can 
be obtained from the equation of the null geodesic, written in terms of 



+ Tlr>Y = , (4) 



where F?. = ^t^^^^-j the Ricci rotation coefficients. Recalling the expression 
of given in equation (2), we finally get 




E 



1 



(5a) 



(56) 



1 



cos - sin 




(5c) 



^2)1/2 



We note that all the information about the spacetime curvature and the flow 
dynamics are contained in the tetrad field and enter the Boltzmann equation 
via the tetrad vectors themselves and their local rates of change which appear 
in the Ricci coefficients. 

Equation (3), together with the set (5), is the more general form of the 
transfer equation and holds for arbitrary flow motions in any given spacetime. 
In the next subsection we will discuss how the existence of spacetime symmetries 
implies that the distribution function is independent on some of the phase-space 
variables, easing the solution of the transfer problem. 



Let us consider the more general spherically-symmetric spacetime, 
described, in spherical coordinates, by the line-element 



Spherical symmetry implies that there exist two constants of the motion, Lz 
and L, which are related to the components of the photon 4-momentum by 
Lz =P3 = sin^ $p^, = r^[(p^)^ sin^ 6 + (p^)^]. These two expressions take 
a very simple form, and lead to a major simpliflcation in the transfer equation, 
if the fluid conflguration and the radiation held are themselves spherically- 
symmetric. In this case the spatial 3-velocity v of the comoving observer. 



b) Transfer in Spherically-symmetric Spacetimes 



ds'^ = gooir, t)dt^ + gii{r, t)dr^ + r^{d9^ + sin^ ddcf)'^) . 
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measured by the stationary observer Sq/^Z—Ooqi is in the radial direction and 
the most convenient choice for the tetrad is 

vV-^oo ^/gii 

el - ( 2_ 0^ 

^"Ia/=^'V^' ' ; (6) 

ei= (0,0,r-\0) 

el = (0,0,0,r-^sin-^^) 

where 7 = (1 — v^)~^^'^. The constants and L may be then expressed in 
terms of the tetrad components as 

Lz = L sin $ sin $ (7a) 
L'' = r^E\l-fj,^). {7b) 

In spherical symmetry, the photon distribution function must be independent 
on both the polar angles </> and 6. Since, from equation (7a), we have ^» = $(^), 
it follows that isotropy in coordinate space implies also that df/d^ — and 
the Boltzmann equation reduces to 

df Q df , df dE df dfi 

In the further hypothesis that the spacetime is stationary, the existence 
of a time-like Killing vector provides a third conserved quantity, po = —E^q, 
which can be used to obtain a simple expression for the photon energy along 
each ray in the LRF 

in the previous expression y = ^^—gao is the specific energy of the fiuid, as 
measured by a static observer at infinity. Clearly the differential operator in 
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the transfer equation is Pfaffian, so it is always possible to solve the Boltzmann 
equation along its characteristic directions, i.e. along the photon trajectories 
in the 7-dimensional phase-space. In the case at hand, these trajectories lie in 
a 4-dimensional hypersurface and can be obtained solving equations (5a) and 
(5b) together with 

dt Q E^o /in \ 

^=p^ = Ey^^,+v). (106) 

Actually, the existence of the two constants of motion L and E^o yields 
analytical expressions for both /x and i?, as functions of r, along each photon 
trajectory: 

^ -Z/^^6^±r(r^ + 6^ffoo)^^^ 
E=-. >^f^!l±lL—^E^ (lib) 



1/2 



y 



where the impact parameter b = L/ E^q has been introduced. Due to spherical 
symmetry, only positive 6's need to be considered, negative values of the angular 
momentum give exactly the specular picture, so in the following b^ will be used 
as a parameter. It can be easily shown that the plus/minus sign in equations 
(11) refers to photons for which ji + v is always positive/negative. This implies, 
see equation (10b), that the radial coordinate is always increasing/decreasing 
along the path and that the condition /x + v = defines the locus of turning 
points for the trajectories. This is just a manifestation of aberration: the 
turning points, in fact, are located where the cosine of the angle between the 
photon and radial directions, measured by the stationary observer, vanishes. 

Specializing to the vacuum Schwarzschild solution, photon trajectories in 
physical space may be divided into three classes (see e.g. Misner, Thorne & 

11 



Wheeler 1973): a) those connecting radial infinity with the event horizon, 
characterized by impact parameters in the range 

< 6^ < 27/4; b) those 
that are trapped in the region 1 < r < 3/2 and c) those for which it is always 
r > 3/2. Trajectories of the latter two types have iP' > 27/4. The limiting value 
= 27/4 corresponds to the circular photon orbit. The plot oi n — iJ,{r; h) and 
of E/Eoo = e{r;b) is shown in figures la and lb for a free-fall velocity law, 
^1 _ 7^-1/2 (,g^^ gggj^ from the the figures, photons starting at the 

horizon can reach infinity with non-zero energy only if they are emitted exactly 
in the radial direction (// = 1) with an infinite energy, while ingoing photons 
that leave infinity with zero angular momentum reach the horizon halving their 
initial energy. At large values of r all rays concentrate at /j, = ±1, as radial 
streaming is approached. Trajectories with an impact parameter equal to the 
critical value 6^ = 27/4 exhibit a saddle point at r = 3/2. 

In the following we will concentrate on the case in which both matter 
and the radiation field are stationary. Under this assumption the distribution 
function depends only on three variables, r, E and and since it is £^ = E{r), 
ji = //(r) (see equations [11a], [lib]), the radial coordinate itself can serve as a 
(non-affine) parameter along the null geodesies. The Boltzmann equation can 
be then integrated in the domain of existence of each photon trajectory. This 
particular choice appears to be convenient for a number of reasons, although it 
poses some numerical problems, as it will be discussed later on. First of all, the 
treatment of boundary conditions is much simpler when the radial coordinate 
is the independent variable and this avoids also the integration of equation 
(10b) along with the transfer equation. Moreover, when scattering is taken into 
account, the source term depends on the integrals of / over angles, which must 
be evaluated at both constant r and E. The knowledge of /(r) avoids the use of 
spline or other interpolation algorithms, which is time-consuming and would be 
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needed in the case of a different parametrization of the photon trajectories. In 
conclusion, at least for what concerns the radiation field, the transfer problem 
can be solved integrating numerically the single differential equation 



# ^ ^9 9 
dr y{n + v)E 



(12) 



for different values of the two parameters h and Eoo- 

At variance with what happens using other methods, like for example 
expansion in PSTF moments, this kind of approach makes a great simplification 
in the mathematical structure of the problem: in fact, the non-grey problem 
can be solved without integration of complicated systems of partial differential 
equations. Moreover, no closure is needed and this formalism gives as result 
the full radial, frequency and angle dependent solution. As we will discuss in 
detail in the next section, is just the knowledge of the angular dependence of 
the distribution function, lost when the moments of / are used as dependent 
variables, that gives the possibility to use the characteristic rays method to 
study the Compton scattering in its more general form; this approach naturally 
preserves the hyperbolic character of the Boltzmann equation. 

In this investigation we focus our attention on the calculation of the 
radiation field and, thus, we restrict our discussion to the case in which velocity, 
density and temperature profiles are fixed a priori, similarly to what was done 
by Mihalas (1980) in the special relativistic case and Schinder & Bludman (1989) 
in the general relativistic, static case. Clearly, the full solution of the radiation 
hydrodynamical problem requires the simultaneous integration of the transfer 
equation together with the Euler, continuity and energy equations that, in turn, 
depend on the gray mean intensity J and on the gray radiative flux H: 



2 ./o J -1 2 jQ j_i 



J ^ - I du I Id/j, = - I du I fv^d/j, 
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2 /"OO /•! -j^ /-oo /•! 

H = - dv I^d/j, = 77 / (iz^ / fu^iidii . 
2 7o ^-1 2 Jo J-i 

The coupled solution of the transfer and gasdynamical equation poses, therefore, 

the same difficulty encountered in the integration of the transfer equation alone 

in presence of scattering. A numerical technique for the solution of the integro- 

differential scattering equation is discussed in section 4a. The same method 

can be applied to the full radiation hydrodynamical problem and an example is 

presented in section 5c. 

c) Boundary Conditions 

Because there is not a one-to-one map between r and ^, equation (12) 
must be integrated twice for each value of 6^, in correspondence with the 
two solutions for ji and E given by equations (11a), (lib). At the same 
time, two difi'erent boundary conditions for the distribution function / must 
be imposed, taking into account that the plus (minus) sign in equations (11) 
corresponds to outgoing (ingoing) trajectories. The boundary condition for 
ingoing characteristics of type a) is prescribed in the standard way: for a non- 
illuminated atmosphere, for example, it is just / = at the outer edge of the 
integration domain. This is also the only condition required to integrate the 
transfer equation along all characteristics of type c), since integration can be 
started at large r with, say, / = and carried out until the turning point 
is reached storing the computed value of /, which is then used as the initial 
condition along the outgoing branch of the trajectory. The remaining rays, 
including characteristics of type b), can be treated much in the same way if 
there exists a region in the flow where the effective depth Tg// becomes larger 
than unity at any frequency and LTE is attained. In this case, in fact, the 
required boundary condition is simply / = B^{T)/E^, By(T) is the Planck 
function at temperature T, at a radius f such that Teff{f) > 1. 
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Although this is the standard case for stellar atmospheres, including 
accretion flows onto compact stars, a different situation may arise when dealing 
with accretion onto black holes: for low values of the accretion rate, for example, 
the flow is optically thin all the way down to the horizon (see e.g. Nobili, TuroUa 
& Zampieri 1991). Now a boundary condition for / must be imposed at r = 1 
for rays starting at the event horizon. Since E goes to infinity there, both the 
distribution function and g must vanish. The product E{ii + v), however, does 
not vanish for all outgoing rays at the horizon, so g = implies also df /dr = 0. 
In order to avoid numerical overflows, integration is started at a radius 
fractionally larger than unity, with the regularity condition df /dr = 0. The two 
rays with 6^ = 27/4 are peculiar since they intersect at r = 3/2 (the saddle 
point) which is also a critical point for equation (12). We still integrate the 
transfer equation along these particular rays taking as a regularity condition 
5r = at r = 3/2. Strictly speaking, this condition is exact only in the case in 
which the effective optical depth is larger than unity at the last photon orbit; in 
other cases there is no physical reason to ask for thermalization and the value 
of / may be undetermined. However, since the radial derivative of / diverges 
at the critical point, we found that, in a finite differences numerical scheme, the 
solution of the differential equation fast relaxes and the final result is probably 
not strongly affected by the value of the distribution function at r = 3/2. 
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3. THE SOURCE FUNCTION 

In the following we deal with an unmagnetized, fully ionized, non-degenerate 
hydrogen gas in which emitters and absorbers are in local thermal equilibrium 
at a temperature T. We consider also the case in which electrons are relativistic 
(T^ 5 X 10^ K), and present a fully general treatment of Compton scattering. 
However, for the sake of simplicity, we focus our attention only on thermal 
emission and absorption together with scattering from free electrons; other 
processes, as pair production and double Compton scattering, that may be 
relevant at such high temperatures, are outside the scope of this paper. In this 
section physical units are used; 7 and r denote the dimensionless photon energy 
and electron temperature, both in units of nieC^; Kp{x) is the modified Bessel 
function of the second kind. 

a) Thermal bremsstrahlung 

The source term for spontaneous emission and absorption, including 
stimulated emission, can be written as 

where r] and x cire the emission and absorption coefl&cients, measured in the 
comoving frame. Because of the assumed equilibrium, Kirchhoff law yields: 

with u = E/h. In the medium we are considering, the dominant true 
emission and absorption processes are electron-proton and electron-electron 
bremsstrahlung; in the following we will indicate as «;// the correspondent total 
opacity. The free-free contribution to the source term is then 
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The photon spectrum from bremsstrahlung is usually described in terms of 
the velocity-averaged Gaunt factor G; in the non relativistic regime tables for 
G have been presented by Karzas & Latter (1961). However, as discussed 
by Gould (1980), contributions to the total energy loss rate due both to 
relativistic corrections in the electron velocity distribution and to electron- 
electron bremsstrahlung are already of order 10%atT~10^K become as large 
as 30 % at T ~ 10^ K. Free-free emissivity from a relativistic thermal plasma 
has been investigated by several authors (see e.g. Alexanian 1968; Quigg 1968; 
Haug 1975; Gould 1980; Stepney & Guilbert 1983; Dermer 1984, Dermer 1986). 
The photon spectrum from e-p emission involves a single quadrature over the 
relative Lorentz factor of the interacting particles 7^ (see e.g. Dermer 1986) 

where daB-H{lilr)/d'^ is the Bethe-Heitler cross section corrected for the 
Elwert factor (see e.g. Heitler 1936) and the number density of 

electrons and protons. The previous expression holds for r <^ nip/ me, so that 
protons can assumed to be at rest in the lab-frame. 

Electron-electron emissivity is more complicated since now both particles 
have the same mass and a quadrupole contribution appears. The standard 
expression involves a five-fold integral of the totally differential cross-section 
(Haug 1975), but, as shown by Dermer (1984, 1986), it can be reduced to 
a triple integral exploiting the covariance of Haug's formula to evaluate the 
cross-section in the CM-frame. The final result is 

2^ /-OO (^^^-^^ r<M 



7.^-1) /•^(^'^)rf7*c^<_,(7*,7r) 



^e-e(7,r)-^^^,^^^/^^y^ '^^'^[2(7.'+l)]^/^io 7* '"^7* 



[2(7r + l)]'/' ^7^ + 7*' 



exp 



277* 
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(see Dermer 1986 for notation). 

The numerical evaluation of both ?7e-p and ?7e-e poses no particular 
problems and has been carried out following Dermer (1986) in the ranges 
2 X 10-2 < T < 10, 2 X 10-2 < hv/KT < 25.12. Numerical results for the 
total Gaunt factor were then fitted with the analytical function (see Stepney & 
Guilbert 1983) 



G = 



{A + Bx)hi{l/x) + C + Dx, x = hu/KT < 2.51 

ax"^ + (3x + 'J + d/x, a; > 2.51, 



deriving, for each r, the set of coefficients A, . . . ,d. The Gaunt factor can 
be then obtained at any value of r and hu/KT by means of a suitable 
interpolation/extrapolation. At temperatures below ~ 10 keV (r^ 0.01), the 
asymptotic limits of Gould (1980) are used for both e-p and e-e emissivity. 

b) Electron scattering 

The second important radiative process we consider is scattering from free 
electrons: we recall that one of the major complications encountered in solving 
the transfer equation comes from its non-local character. In fact, even limiting 
to the coherent and isotropic case, the source term is 



^ = Ql^es {ju - /) , (15) 



where Kes is the Thomson opacity and 



1 

> = 2 y ^ / ^) ' 



is the zero-th moment of the distribution function. Allowing for the more 
realistic case of Thomson scattering, the correspondent cross section has a 
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monopole plus a quadrupole angular dependence (see e.g. Chandrasekhar 1960) 
yielding 

(16) 



Qes 

E 



QKe 



where 



1 

2 J ^ ^'^^^ ' 



The Thomson limit can be assumed to correctly describe electron scattering 
when the energy exchange in a single collision can be safely ignored. On the 
other hand, in high temperature regions non-conservative effects and quantum 
corrections play a fundamental role in shaping the emergent spectrum. The 
derivation of the general expression for the Compton source term can be found 
e.g. in Pomraning (1973) and is briefly outlined below, mainly to introduce 
some basic ideas which will be used later on when the numerical scheme is 
discussed. With reference to a single scattering, n denotes the incident photon 
direction and ^ = ft - ft' , where primed quantities refer to the scattered photon. 
For an incident photon energy 7 and an electron velocity the Klein-Nishina 
formula gives the probability of scattering into the energy 7' and the direction 



n 



(7(7 



/ -> 

7 ,n 



n ,Ve) 



1 + 



(1-0 
X^DD' 



+ 



(1-0 

X^DD' 



:i7) 



X 8 



D 



D' 



i-i + x--x— 



7 



7 



where 



D = 1 — n ■ Ve/c, D' = 1 — n' ■ Ve/c, 



X 



1-4 



-1/2 



ro is the classical electron radius and S is the Dirac delta function. Integration 

over the relativistic maxwellian distribution 

A^exp(-A/T) 



fe{Ve) = 



4:7rTC^K2 (l/r) 
19 



(18) 



gives the Compton Scattering Kernel (CSK) 



dv, 



fe {Ve 

^ A 



1 + 



X^DD' 



+ 



(19) 



Here the CSK is normahzed with respect to KesQ, which is reciprocal of the 
Thomson mean free path; the inverse probability, related to the scattering 
emissivity, can be obtained from the detailed balance condition 

(7 ^ 1', ^, t) 7^ exp (-7/t) = ct (7' ^ 7, ^, r) 7'^ exp {-i /t) . (20) 

Further integrations over all outgoing photon directions and energies 
provide the source term appearing in the Boltzmann equation 
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Inserting equation (20) into equation (21), the latter can be written in the more 
compact form 
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where /' = /(r, n', 7') and 

aoo= / / dnV(7^7',e,T) 
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(22) 



(23) 



is the zero-th moment of the CSK (Shestakov, Kershaw & Prasad 1988). In 
the previous expressions, non-linear terms account for stimulated scattering. 
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The general task of computing the moments of the CSK was undertaken by 
Shestakov et al. (1988). They have shown that by performing the integration 
over 7' first and exploiting the 5-function, the expression of the zero-th moment, 
which is originally a fivefold integral, can be reduced, after a lot of non-trivial 
algebra, to a single quadrature 



where: 



y^o {y) =^ 
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(246) 
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The full evaluation of the Compton source term involves a number of 
very complicated six dimensional integrals of the distribution function weighted 
by the CSK for each value of 7, r, ji. Because only discrete values of the 
distribution function will be available, all the six quadratures should be, in 
principle, evaluated numerically at each grid-point and this would make the 
integration of the transfer equation prohibitively time-consuming. However, 
as discuss by Kershaw, Prasad &; Season (1986), two of the three integrals in 
the CSK become analytical if a particular polar axis for projecting the electron 
velocity is chosen. Moreover, Kershaw (1987) presented an efficient method 
for calculating the single integral of the CSK over 7' or ^ and the double 
integral over both these variables. A detailed discussion of our algorithm for the 
evaluation of the first addendum in the Compton source term, that is essentially 
a re-adaptation of Kershaw's method, is presented later on. 

Although the treatment we have just described is the more general to 
handle Comptonization and proved to be reasonably fast, it remains very time- 
consuming, so it is useful to have approximated expressions of gc that can be 
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used in some regimes. As it is well known, the complicated nature of the CSK 
has led many authors to model the Boltzmann equation by a diffusion equation 
in the frequency space. This approach, the Fokker-Planck approximation, was 
firstly used by Kompaneets (1956) in the limit of small 7 and r. Relativistic 
corrections to the Kompaneets equation can be included modifying the diffusion 
coefficient, and a number of efforts were devoted to extend its original form 
(Eraser 1966, as quoted in Pomraning 1973, Cooper 1971). More recently, 
Prasad et al. (1988) derived an exact analytical expression for the diffusion 
coefficient that holds for arbitrary values of 7 and r, in the assumption of 
a nearly isotropic radiation field. The main simplification introduced by the 
Fokker-Planck approximation is that the integral operator in the transfer 
equation is replaced by an infinite order differential operator that, for small 
values of 7 and r, truncates at a finite order. The method, originally developed 
for the non-relativistic transfer equation, is based on an expansion of the specific 
intensity in a Taylor series about ly' = v. At the first order in 7 and r, Eraser's 
result is 



gcE^ = - KesQ (1 - 27) / + KesQ [ dVt' f Pn (0 SnI 



(25) 



where is the Legendre polynomial of order n and Sn {n = 0, . . . , 3) are 
second order differential operators (see Pomraning 1973). Using the standard 
relation 



e = + Vl-/^'Vl-/^^cos ($ - $') , 
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the previous expression can be cast into the form 



^ =KesQ [Ai + IJiA2 + (1 - I?) ^3 + (3 - ^1?) AA - KesQf {1 - 2^+ 

hi 

T [^5 - ix^A^ + IX (3//2 _ 5) ^7 + ^ (3 _ 5^2^ As] } . 



The quantities Ai^ containing the first four moments of /' and their first 
and second frequency derivatives, are reported in Appendix A. This is the 
expression of gc needed in the general relativistic transfer equation in Eraser's 
approximation. We stress that up to now no assumptions have been made about 
the angular dependence in the energy exchange terms. A further simplification 
can be introduced if all terms, but /, in equation (26) are assumed to be 
isotropic and are replaced with their zero-th moments. The Compton source 
term becomes then 



is the mean intensity. The approximated expressions (26) and (27) are to 
be preferred whenever a non-relativistic plasma is considered, since their 
evaluation is much faster than that of the general source term given by equation 
(22). Moreover, equation (27) contains far fewer terms than (26), and has the 
great advantage that all the angular dependence is contained in /. 

All forms of the Compton source term based on the Fokker-Planck 
approximation contain both first and second frequency-derivatives of the 
moments of the distribution function. As noted by Nobili, TuroUa & Zampieri 



(26) 




(27) 




where 




23 



(1993), in connection with the system of the first two PSTF moment equations, 
Compton terms act as singular perturbations, changing the mathematical 
character of the differential operator that becomes elliptic. As we discuss in 
detail later on, our numerical code is based on an iterative scheme in which 
integral terms, together with their derivatives, are treated as forcing terms, the 
only full-fledged differential operator being the one contained in the Boltzmann 
equation. On the other hand, the characteristic ray method provides the angular 
and frequency dependence for / that allows to write the Compton source term 
in its original form without resorting to the Fokker-Planck approximation. In 
this case the problem of radiative transfer with comptonization can be solved 
exactly in any range of energies and optical depths, and the hyperbolic character 
of the Boltzmann equation is naturally preserved. 

4. THE NUMERICAL METHOD 

In this section we describe in some detail the numerical scheme we have 
developed for solving the transfer problem. The more general case, which 
corresponds to spherical flows in a Schwarzschild spacetime, is discussed in 
subsection a); in subsection b) a simplifled version of the code, for the solution of 
the full radiation hydrodynamical problem in static, plane-parallel atmospheres 
is presented; flnally subsection c) is devoted to the numerical evaluation of the 
Compton source term. 

a) The spherical case 

As it is well known, in a scattering medium, the transfer equation is an 
integro-differential equation, while it has a simple structure when only true 
emission-absorption is included; in particular, it reduces to an ODE when 
written in its characteristic form. This suggests that its solution can be found 
using an iterative method in which the starting point is just the solution of the 
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transfer problem with only free-free processes taken into account. Following 
this idea, equation (12) has been integrated numerically, with the boundary 
conditions previously discussed, for a given set of values of the parameters 
and E(x> and with the source term g = gff. This provides the zero-th order 
approximation, f'^^\r,iJi,E) of the distribution function, that can be used to 
evaluate the scattering integrals appearing in g^g or gc- In the second step, we 
use the full expression for g to obtain the first order approximation f^^\r, E). 
This is the solution of the transfer equation written in the form 

^^^^ = f + «[/<°'l-«/<V". 

All the expressions of the scattering source term discussed in the previous 
section can be cast, and have been presented, in this form. In equations (15), 
(16), (26) and (27) (3 can be immediately identified with the coefficient of /; in 
equation (22) a is the integral term. The scheme is iterated until convergence is 
reached, improving at each iteration the functionals a and P making use of the 
distribution function computed in the previous step. As a convergence test, we 
compared each element of the matrix j^, with its value relative to the previous 
iteration and stored the maximum relative correction. Cauchy criterion has 
been applied to verify the convergence of the succession of such corrections. 

Equation (12) has been integrated using a a finite differences method 
originally developed by Nobili & TuroUa (1988), in which the algebraic system 
is iteratively solved using the Henyey technique for matrix manipulation. The 
entire radial domain [rin,rend] is divided by M points; rays of type a) are 
integrated using this grid. For trajectories which exhibit a turning point, the 
transfer equation is solved on the same mesh, picking up the subset of grid points 
which cover their region of existence. Although, as we already mentioned, the 
choice of r as the parameter along the geodesies has a number of advantages, it 
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results in a divergent derivative of / at = —v. While, for those branches which 
approach the turning point, this introduces some errors at most in the last few 
points, trajectories moving away from the turning points may be systematically 
affected by an inaccurate determination of their boundary condition. However, 
it should be taken into account that when the optical depth at the turning point 
is either large or very small, / tends to Bj^/hcE^ or remains vanishingly small, 
independently on the boundary condition for equation (12). Numerical errors, if 
any, are, then, restricted to rays inverting in regions of moderate optical depth. 

The choice of the 6-grid strongly constraints the final angular resolution 
of /, and requires special care. Let us first assume that a black hole is the 
central source; in this case the interval < 6^ < 27/4 corresponds to ingoing 
and outgoing trajectories of class a). For these two subclasses, we fix Ni and 
A^2 values of the impact parameter in such a way to produce an equally-spaced 
(U-grid, in the range [—1,1], at the critical point r = 3/2. To discretize the 
range 6^ > 27/4, we exploit the one-to-one correspondence between and the 
position of the turning points, r = Tn, 

6^ = ^ . (28) 

We fix and A''4 values of rn, the first at r < 3/2 and the latter at r > 3/2; the 
Tn's are just the radii of the spherical shells tangent to the orbits of types b) and 
c). In such a way, the total number N of //-points in the interval — 1 < /x < 1 is 
r-dependent, and it is bounded by iVi + iV2 + 1 < < iVi + + SiVs - 1 for 
r < 3/2 and N1 + N2 + I <N< N1 + N2 + 2N4-I ior r > 3/2. A better angular 
resolution in all the radial domain can be obtained increasing the number of 
photon trajectories. In the case the central source is a star of radius r*, the 
/i-gridding works in a very similar way, but the values of b in the range 



r* - 1 
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now produce an equally-spaced //-mesh at the star radius, the A^4 points refers 
to r > r* while no trajectory of type b) is present. We have found more 
convenient to derive the values of the impact parameter starting from the radial 
coordinate of the turning points, and not vice versa, since in this way the radial 
extent of the photon trajectories, and hence the integration range of equation 
(12), is specified without solving the cubic equation (28). 

Once the rays are fixed, equation (12) must be integrated for different values 
of the parameter Eoo along each trajectory. The range of £^oo should be chosen 
in such a way that, at each value of r, we can compute the distribution function 
in an interval of the local energy, [Emin, Emax], large enough to cover the 
interesting portion of the spectrum. The parameter range [(-E'oo)min, (-E'oo)maa;] 
must be larger than [Emin, E^ax] at any given radius, since both gravity 
and dynamics act in changing the photon energy along the geodesies. For 
1^ < fend^ in fact, the energy interval [Emim E^ax] is actually influenced by 
some characteristic rays starting at Tend with E^o in the range 

{Eoo)min = [y (1 + ^)]r^,„ Emin < -^oo < (1 - ^)]r^i„ E^ax = {Eoo)max ■ 

In the numerical calculations we have used the dimensionless energy x = 
E/KT^, where T* is a suitable normalization temperature. For later appli- 
cations, we found more convenient to divide the storage window [xmini^max] 
by means of L points equally-spaced in Inx; the same grid is maintained at 
all radii and / is stored at these points as a function of the local dimensionless 
energy using an interpolation. In the two remaining ranges [{Xfx^)miniXmin\ 
and [xmaxi{xoo)max\i 2P valucs of a^oo has been specified. For these values 
of the energy at infinity, the transfer equation has been integrated only along 
the trajectories of those photons that, at some r, have a local energy within 
our storage window. Loading the matrix /(r^, Ek) is particularly convenient 
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since it allows a more direct calculation of the scattering integrals, that are 
evaluated at both constant r and E. All angular integrals can be obtained 
simply performing a weighted sum of / over the //-index without any additional 
scanning of the array or extra interpolations. This has, also, the advantage that 
the we are free to choose the most suitable numerical scheme to integrate over 
energies since the rearrangement of the energy points at each radius follows 
automatically. The numerical evaluation of the frequency-dependent moments 
of / has, however, to be carried out with some care. In particular, when 
the optical depth drops below unity and radial streaming is approached, the 
integration over fi becomes troublesome and we found more convenient to 
perform the quadrature over b^, using equation (11a). Since the same change of 
variable works well near the horizon, where outgoing rays concentrate towards 
= 1, it has been used in all the radial range. However, because of the 
divergence of d/i/db^ where /i = —v, in a small region around this point the 
original //-integration was performed at each value of r. 

b) The Static, plane-parallel case 

The numerical scheme we have just presented allows the solution of the 
transfer equation along the geodesic rays in the more general case, when 
gravity, dynamics and sphericity are all accounted for. In many astrophysical 
applications, however, transfer of radiation through a static, geometrically- 
thin atmosphere is of interest, like, for example, when studying reprocessing of 
thermal radiation in the atmosphere of X-ray bursting neutron stars. In all 
these plane-parallel approach to the solution of the transfer problem is 

fully justified since the atmospheric scale height is much less than the star 
radius, although the effects of the strong gravitational field must still be 
considered. The assumption of hydrostatic equilibrium introduces a major 
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simplification in the treatment of radiative transfer because advection and 
aberration are no more present. For a vanishing velocity field, equation (9) 
reduces to E = £^oo/-\/~^00) implying that the value of the local energy 
at a given radius is the same along all rays. This is just another way of 
stating the existence of Thome's (1981) Universal Red-shift Function. The 
rays are now symmetrical with respect to the fj, = line. A further, drastic, 
simplification follows if it can be assumed that the radial coordinate is constant 
in the atmosphere and equal to the star radius. This is commonly done in 
non-relativistic transfer theory, replacing the height above the base of the 
atmosphere with the optical depth. The rays are just straight lines, /i = const, 
while the photon energy seen at infinity is simply the energy at any depth red- 
shifted by the constant factor (1 — 1/r*)^/^. In the light of these considerations, 
we have developed and tested a reduced version of our code which uses the 
scattering depth as the independent variable. The angular mesh is obtained 
specifying directly the values of the energy points at which / is computed 
coincide with the energy grid, which is the same at all depths. The calculation 
proceeds exactly in the same way as in a non-relativistic problem and the 
spectrum at infinity is simply obtained by applying the gravitational red-shift 
factor to the spectrum emerging at the top of the atmosphere. 

An application to isolated neutron stars accreting at low rates is presented 
in section 5b. In this problem electrons are far from being relativistic 
so Comptonization can be safely treated in the diffusion approximation 
using expression (27) for the scattering source term. The much shorter 
computational time allowed us to solve also the thermal and pressure structure 
of the atmosphere, coupling the hydrostatic balance and the radiative energy 
equilibrium to the transfer equation. The hydro equations are solved iteratively, 
exploiting the scheme for the computation of the scattering integral we have 
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already discussed. Pressure and temperature profiles are computed at each 
iteration step, once the frequency-integrated moments have been obtained. 

c) Numerical evaluation of the Compton source term 

As discussed in section 3, the Compton source term, in the form (22), is the 
sum of two contributions. The second addendum, which requires the calculation 
of the zero-th moment of the CSK, aoo, poses no problems since it involves a 
single quadrature of an analytical function. As proposed by Shestakov et al. 
(1988), upon the change of variable 

(Too can be efficiently evaluated using a Gauss-Hermite quadrature. We have 
tested that six points give an accuracy better than 3 parts in 1000, sufficient 
for our purposes. 

We are left, then, with the problem of finding a fast algorithm for the 
numerical calculation of the multiple integral 

/ dl]V(7^7',e,T)/'. (29) 

First of all, we note that the scattering probability may become strongly peaked; 
in the Thomson limit, for example, the CSK tends toward a 5-function at 
7 = 7'. In all regimes in which the integrand is fastly- varying particular care 
must be used to account for delicate cancellations between opposite terms. We 
start considering the CSK itself. As discussed by Kershaw, Prasad and Beason 
(1986), the complicated three-dimensional integral in the electron velocity space 
can be reduced to a single integral when the solid angle element is defined with 
respect to a particular polar axis. In fact, taking the polar axis in the direction 
of the photon momentum transfer s = (7'n' — 7n) /q, q = •\/7^ + 7'^ — 277'^, 
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and using the Dirac 5-function to integrate over the polar angle, the integration 
over the azimuthal angle becomes analytical. The final form of the CSK is then 

^ ^' ^' = 327.rl. (1/r) ^'^^^^^ { ^ 
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where the Lorentz factor A is now the integration variable, o;^ = (1 + ^) / (1 — ^) 
and 
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As stressed by Kershaw et al., the main features of the scattering 
probability are contained in the exp{—X+/T)/q term: everything is smoothly 
varying with respect to this quantity, in particular with respect to the 
exponential. Kershaw et al. proposed two methods for the numerical 
evaluation of the A-integral in equation (30); in both cases the CSK is 
reduced to an approximate analytical expression. Here we adopt their fastest, 
although less accurate, algorithm which is based on a suitable division of 
the integration domain into subintervals where the exponential is replaced 
by a linear interpolation. To avoid delicate cancellations when r ^ 0, a 
Taylor expansion of the inner expression in curly brackets is used to obtain 
an asymptotic series in terms of Legendre polynomials for the integral; only 
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terms up to second order are retained. Using this method the evaluation of 
the CSK becomes analytical with an accuracy of about 3 parts in a thousand 
in all parameter ranges. The CPU time for a single evaluation is typically few 
microseconds on an alpha DEC-3000. 

The algorithm we adopt for computing integrals involving the CSK follows 
the original method presented by Kershaw (1987) for evaluating the total 
scattering cross-section and it is based on the fact that A+ has a minimum 
in both 7' and ^. The most important contribution to the CSK comes, in fact, 
from regions near this minimum; everywhere else the scattering probability 
goes to zero exponentially fast with an e-folding length that is simply r in A+. 
Having these considerations in mind, the double angular integral in expression 
(29) can be written, taking ^ and as the polar and azimuthal angles, as 

/ dn'a (7 ^ 7', r) /' = d^a (7 ^ 7', ^, r) * (r, fx' , 7') . 

J4n J-1 Jo 

For each value of r, /i, 7', ^ the azimuthal integral is evaluated using a Lobatto 
quadrature. The values of the distribution function at 

where 0^ are the Lobatto abscissae, are obtained from a linear interpolation. 
Once this is done, for each value of 7, 7', r, the integration over all polar 
directions is carried out picking up the ^ range, within the interval |^| < 1, 
that provides a non-negligible contribution to the scattering probability: as we 
anticipated, this is the region around the minimum of A_|-. For fixed 7, 7' and 
T, the e-folding lengths in A+, nr, immediately provide the e-folding lengths 
in ^, Denoting, in fact, with = 1 — I7 — 7'|/(77') the value of ^ where 
A+ is minimum, ^„ is the root of the equation 

A+mi + nr = A+ (7, 7', Cn) , 
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where A+^i = A+ (7, 7', Cm) and 

Cm = max (min (Cm, -1) , 1) • 

Within each e-folding interval we use a 4-points Lobatto quadrature and the 
number of intervals is fixed by the request that either the fractional contribution 
of the last e-folding is less than the desired accuracy (5 x 10~^ in the present 
case) or the boundary C = ±1 is reached. At 7 = 7' and C = 1 the CSK has an 
integrable (~ \/l — C) singularity, that can be easily eliminated with the change 
of variable 77 = \/l — C- 

The integration in energy is carried out in a similar way. Since the most 
important contribution to the inner integral (over C in our scheme) comes from 
regions where A-|- is near A+rni, the larger contribution to the outer integral 
(over 7') is provided by regions where A+^i itself is minimum. Clearly, the 
lowest values of A+mi correspond to Cm = Cm, i-e. to Cm > —1; the inequality 
Cm < 1 is always satisfied. We distinguish two cases: for 7' < 7 the previous 
condition is verified in the interval 

— - — < 7' < 7 
1 + 27 - - 

that we call region A, while for 7' > 7 it holds in two different domains, that 
we call in general region B, depending on the value of 7: 

l<l'<Y^ if 7<l/2 
7<7'<oo if 7>l/2. 

Since in region A A+mi = 1, the search for the e-folding lengths is not required 
and integration is straightforward. This is not the case in region B, where 
A+mi = 1 + 7^ — 7. Now, although its minimum value is still A+m2 = 1, A+mi 
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is not a constant. The corresponding e-folding lengths 7^ are to be derived 
solving the equation 

A+m2 + nr = A+ (7, 7^,, ^m) , (36) 

which reduces to 

l + nT=l + 7;-7 (37) 

and gives simply 7^ = 7 + nr. In region B integration over 7' is carried out 
using the same procedure introduced for the ^-quadrature. To complete our 
discussion, we need to consider the two intervals 

region C, and, if 7 < 1/2, 

region D. In both cases X+mi = X+ and its minimum is reached at 

or 

A+m2 = A+ml|-|,/=^/(i_2^) 

in regions C and D respectively. The corresponding e-folding lengths are 
obtained from equation (36). Lobatto rule is used everywhere and stepping 
is terminated when its fractional contribution becomes less than 5 x 10~^. 
The most convenient number of Lobatto points depends on the typical relative 
values of the photon energy and gas temperature. In fact, for different values 
of 7 and r, A+(7') can be either strongly peaked or very broad near its 
minimum. Optimization requires some numerical experimenting, looking for the 
best agreement between the direct evaluation of the CSK double integral and uoo 
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computed using equations (24). For the test model presented in subsection 5b, 
we used either a six or a ten points quadrature. Accordance between the values 
of (Too obtained using the two methods is better than 3-4%, with the larger 
errors in the external region where the radiation temperature (mean photon 
energy) is very far away from the gas temperature. On the other hand, where 
the Compton parameter Yc (see e.g. Rybicki & Lightman 1979) is greater than 
unity, accuracy is better than 7 parts in a thousand. We finally note that the 
choice of a gaussian-type quadrature is motivated, basically, by the fact that we 
need to perform integrals of the CSK times /. The distribution function must 
be interpolated to obtain its values at the integration points. Clearly, gaussian- 
type quadratures with a fixed number of abscissae are much faster, although less 
accurate, than step-adaptive schemes, as the Simpson rule originally used by 
Kershaw. Computational feasibility is also the reason for which we decided to 
evaluate the integral (29) using the values of / relative to the previous iteration. 
Clearly, it is possible to rewrite expression (29) as 



since / does not depend on 7'. Now / is just the dependent variable of the 
transfer equation at any iterative step. The drawback is that the computing 
times is about doubled, because of the two multiple integrals. The CPU time for 
a single evaluation of expression (29) is typically ~ 0.1 s and, in an production 
run, a ~ 2 X 10^ evaluation are required, implying a total time of about 6 hr. 
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5. APPLICATIONS 

As we stressed several times, accretion flows onto compact objects, and black 
holes in particular, provide an ideal arena for applications of relativistic 
radiation transfer in differentially-moving media. The accreting matter reaches, 
in fact, not only r ~ 1 with v ~ 1, but, often, temperatures high enough 
(T^ 10^ K) to make a full treatment of Comptonization necessary. For this 
reason we decided to present in this section the numerical solutions of the 
transfer problem relative to different accretion regimes onto black holes and 
neutron stars: "cold" and "hot" black hole accretion is considered in subsections 
a) and b), respectively, and subsection c) deals with "cold", static atmospheres 
around neutron stars. The full radiation hydrodynamical problem in solved only 
in the latter case, while in the first two examples the flow hydrodynamics is kept 
fixed. Since our present goal is to test the capabilities of our method, no attempt 
has been made to explore the models parameter space: we just present results 
for a single model which we judge useful in illustrating the main features of our 
integration scheme. For model c) a direct comparison with the results obtained 
by Zampieri et al. (1995) with the moment expansion has been made, showing 
a good agreement. No previous solutions for black hole accretion spectra are 
available, at least for models which contain an optically thick core. Our attempt 
to cross-check results presented in section 5a integrating the moment equations 
were hindered by severe numerical problems which arise when the fiow is not 
effectively thick at the horizon at all frequencies. The moment method can not, 
also, be used to compute radiative transfer in "hot" models, where Compton 
scattering must be treated outside the Fokker-Planck approximation. 

a) Accretion onto black holes: low-luminosity solutions 
Spherical, stationary accretion onto a Schwarzschild black hole has been 
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throughly investigated in the past and we refer to the paper by Nobih, TuroUa 
and Zampieri (1991, NTZ in the following) for all details. A distinctive feature 
of black hole accretion is that, for the same value of the accretion rate m 
which is the only free parameter, two solutions may exist with very different 
properties: a "cold", low-luminosity and a "hot", high-luminosity one. In 
both solutions the sonic point is so far away that we can safely assume that 
matter is free-falling with = r~^/^ in our region of interest. Here we refer to 
models with high accretion rates, m> 1 in Eddington units. In this regime 
low-luminosity solutions start to develop an inner region optically thick to 
both free-free and scattering and show negligible Comptonization; consequently, 
electron scattering can be treated in the Thomson limit, using expression (16). 
Under these conditions we expect, however, bulk motion Comptonization in the 
converging flow (Blandford & Payne 1981; Payne & Blandford 1981; Nobili, 
TuroUa & Zampieri 1993) to act efficiently at high frequencies where true 
absorption is very low. 

For our first test we consider NTZ solution characterized by m = 0.71, 
corresponding to a density at the horizon qh = 10~^gcm~^. The gas 
temperature is in the range 2 x 10^ 5 x 10^ K, so we have chosen a 

normalization temperature InT* = 11. The dimensionless energy window is 
a^min = 0.1 < a; < Xmax = 40, Corresponding to the range 0.5-206 eV; here 
L = 30 points have been used. Outside this range, P = 10 energies has been 
fixed in each of the two additional intervals of Xoo we need, as discussed in 
section 3a. Since the effective optical depth at our larger energies is everywhere 
< 1, we solved the transfer problem for 10~^ < log r < 5, imposing the boundary 
condition df /dr = along trajectories starting at ri^. The radial domain has 
been divided by M = 250 points; the grid is not uniformly spaced and points are 
tighter around r = 3/2. To obtain a good angular resolution, 90 trajectories 
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have been followed at each energy, AT^ = A^2 = 20, A^s = 10, = 40. In 
such a way, the number of // points, which is minimum at r — 3/2, is always 
greater than 21. At each value of x within the storage window the scattering 
source term was calculated using a linear interpolation for both the matrices 
j,^ and kiy] outside this window an extrapolation has been used. Figures 2a 
and 2c show the mean intensity J^, and radiation pressure at different 
energies, together with the Planck function at T{rin); each curve corresponds 
to a different value of the radial coordinate. The effective frequency-dependent 
optical depth goes from 3 x 10^ to 10""^ for the lowest frequency, while high 
energy photons stream freely at all radii. The low energy portion of the spectral 
distribution is a superposition of thermal bremsstrahlung emission at different 
temperatures while bulk motion Comptonization produces a power law high- 
energy tail. The calculated spectral index, a = —2.9, is due to unsaturated bulk 
motion comptonization, being the scattering optical depth ~ 0.7 at the horizon. 
The theoretical value derived by Payne & Blandford, in the limit t^s ^ 1, 
is Q! = — 2 for a free-fall velocity. At large radii, where radial streaming is 
approached, all moments fall off as r~^; the asymptotic radial gradient we have 
found is —1.99. While the evaluation of even moments does not pose particular 
problems, in the inner regions, where the radiation field is nearly isotropic, a 
direct numerical quadrature for computing odd moments becomes troublesome 
because of the delicate cancellations between contributions of opposite sign. To 
avoid this problem, the monochromatic flux, presented in figure 2b, has been 
replaced with its analytical expression in the diffusion approximation every time 
it is Teff > 10. A typical production run required 10-11 iterations to converge 
with a fractional accuracy better than 10~^, with a total CPU time of about 
20 minutes on an alpha DEC-3000. 

b) Accretion onto black holes: high-luminosity solutions 
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In "hot" solutions temperature is much higher, typically ~ 10 ° K near 
to the horizon. As a consequence, free-free absorption is much lower than 
in "cold" models, even for larger accretion rates. Along the high-luminosity 
branch, thermal Comptonization is the dominant radiative process and it must 
be treated in its more general form, using expression (22). Here we consider 
the "hot" solution of NTZ with m = 71, qh = lO'^gcm-^. The flow is 
now effectively thin at all frequencies, although an inner core optically thick 
to scattering is present. The gas temperature in this model is in the range 
10^ K < T < 10^^ K, so we have chosen InT* = 21 and x^ain = 0.008. Now 
the energy window is 0.9-4500 keV and L = 35 points have been used. Since 
the evaluation of the CSK integrals is very time-consuming, both angular and 
radial resolution has been reduced with respect to the previous model: A^i = 
N2 = = 10, A^4 = 30 and M = 110 in the same radial domain. In this model 
convergence has been reached with a fractional accuracy better than 0.02, and 
the calculated radial gradient at infinity is —2.04. The resulting mean intensity 
is presented in figure 3. In high temperature models, the mean intensity is 
always less than B^, but despite the accreting gas radiates less efficiently than 
in low-luminosity optically thick solutions, the efficiency of accretion process 
is higher, due to the fact that the matter temperature is now higher in the 
whole photospheric region. Since the emergent spectrum is peaked at about 40 
keV, these solutions, if stable (see NTZ and Zampieri, Miller & TuroUa 1995), 
seem to provide a natural way to produce hard X-ray radiation with reasonable 
efficiency out of spherical accretion onto black holes. 

c) Static, plane-parallel atmospheres around neutron stars 

Our last application refers to a static, "cold" atmosphere around an 
accreting unmagnetized neutron star. Emitted spectra were firstly derived by 
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Zel'dovich & Shakura (1969) and, in more detail, by Alme & Wilson (1973) 
in the high luminosity range (Z^ 10~^ in Eddington units). Solutions for 
10~^ < I < 10~^ has been recently presented by Zampieri et al. (1995) in 
connection with isolated neutron stars accreting the interstellar medium. As 
they have shown, the emitted spectrum exhibits an overall hardening with 
respect to the blackbody at the neutron star effective temperature with a 
hardening ratio, typically ~ 1.5 — 3, which increases with decreasing luminosity. 
The most important physical processes are free-free emission and absorption; 
Compton cooling plays a role only in increasing the temperature in the outer 
atmospheric layers, where the low energy tail of the emitted spectrum is 
created. Since typical temperatures are very low, Comptonization is treated 
by means of the approximated expression (27), which is much faster. The run 
of thermodynamical variables is obtained solving the hydrostatic and energy 
balance (see Zampieri et al.) together with the transfer equation, using an 
iterative scheme. However, in this problem the thermal balance is very delicate 
and the zone where photons of different energies thermalize strongly depends 
on integrated quantities, J and the absorption mean kq. Numerical integration 
proved more stable if J and H are derived as solution of the first two gray 
moment equations. The same approach was used by Zampieri et al., with the 
difference that in our scheme the gray moment equations are solved exactly, 
computing from the specific intensity the Eddington factors K/ J ai each depth 
and J/H at r = 0. Here we have recomputed the model with I = 10~^, using 
a spectral window Xmin = 0.1, Xmax = 10 centered around a normalization 
temperature logT* = 6.6; L — 30 frequency points have been used. The angular 
resolution is provided by Ni = N2 = 15 trajectories and the transfer equation 
is solved in the range —8 < log(T) < 0.9 using M = 100 grid points. The 
resulting mean intensity is plotted in figure 4. Figures 5 and 6 show the emergent 
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spectrum and the temperature profile together with the results of Zampieri et al. 
(dashed lines). As can be seen, the agreement both in the spectral shape and the 
temperature profile is very good, showing that the approximated solution of the 
transfer equation with two moments is rather accurate. To obtain this model, 
with a fractional accuracy better than 2 x 10~^, 14 iterations were required, 
with a total CPU time of about 3 minutes on an alpha DEC-3000. Agreement 
between the gray mean intensity, derived as the double integral of /, and the 
solution of the second gray moment equation is always better than few parts in 
thousand. 

6. CONCLUSIONS 

In this paper we have presented a characteristic method for the solution 
of the general relativistic transfer equation. If the spacetime admits some 
symmetries, the formalism can be simplified; in particular, in presence of three 
Killing vectors, two of the three equations for the characteristic rays become 
analytical. In addition, using the radial coordinate as the parameter along 
the null geodesies, the exact solution of the transfer problem can be obtained 
solving a single ordinary differential equation along a bi-parametric family of 
characteristic trajectories. A numerical technique, based on an iterative scheme, 
has been developed and tested either for the calculation of the radiation field 
in a fixed background or for the solution of the full radiation hydrodynamical 
problem in spherical and plane-parallel geometry. Particular care has been 
devoted to the evaluation of the source term, taking into account radiative 
processes which are believed to be of importance in astrophysical accreting 
plasmas: electron-electron and electron-proton bremsstrahlung, Thomson and 
Compton scattering. 

Radiative effects due to magnetic fields and pair production-annihilation 
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were not considered in this work. However we stress that the method we 
have presented is completely general and additional radiative processes may be 
easily included. Source terms not involving integrals of the photon distribution 
function can be simply accounted for when the corresponding emissivity and 
opacity coefficients are provided. On the other hand, our iterative scheme allows 
for the solution of integro-differential equations and can be used to include 
also different integral source terms as, for instance, those ones related to pair 
production or bound-bound emission. Actually, a self-consistent treatment 
of pair production entails the solution of the full radiation hydrodynamical 
problem, with the addition of the pair balance equation and was left out on 
purpose in our discussion of "hot" accretion solutions which are obtained at 
fixed hydrodynamics. 

In the test models we discussed magnetic fields and pair production are not 
expected to play a relevant role at least for low luminosity solutions. In the case 
of accretion onto neutron stars it can be easily shown, in fact, that, for typical 
temperatures and densities in the photospheric region, the cyclotron emission 
is lower than the free-free emission if 10^ G (see e.g. Schmid-Burgk, 1978). 
On the other hand, a relic magnetic field of this order is just what is expected in 
isolated neutron stars which evolved beyond the pulsar phase; our models can 
be then assumed to describe correctly the emitted spectrum from old neutron 
stars accreting the interstellar medium. As far as low-luminosity accretion 
onto black holes is concerned, the limiting value is a factor 10~^ smaller, but 
it still exceeds the maximum strength of the tangled S-field derived assuming 
equipartition between magnetic and thermal energy densities. However, in the 
inner regions of high-luminosity models electrons become relativistic and both 
pair processes and synchrotron emission start to play a role. 
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APPENDIX A 

Here we present the expressions of the Ai, i = 1, . . . , 8, terms appearing into 
equation (26). Their derivation starts from Eraser's result, equation (25), and 
makes use of relation 



At first order in 7 and r, only the first four moments of the distribution 
function 



appear in the Compton source term; they will be termed hj^^ kj^ and l^,^ 
using the standard notation. We introduce also the correspondent moments of 
the specific intensity /: 



To make the larger number of terms dimensionless, we introduce double 
logarithmic frequency derivatives for all even moments; terms containing odd 



i = ^in' + Vl-A^'Vl-/^^ cos ($ - $') . 
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moments are written as semi-logarithmic derivatives since odd moments may 
become vanishingly small in regions where the effective optical depth is very 
large. The final result is: 
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FIGURE CAPTIONS 

Figure la. The run of the cosine of the angle between the photon momentum 
and the radial direction, as measured by a free-falling observer, along the 
characteristic rays. Different curves correspond to different values of the 
impact parameter b. 

Figure lb. Same as in figure la for the photon energy normalized with respect 
to Eoo- 

Figiue 2a. Monochromatic mean intensity at different radii (full lines), 
together with the blackbody function at T(rj„) (dashed line), for "cold" 
accretion onto a black hole with m = 0.71. 

Figure 2b. Same as in figure 2a for the monochromatic flux. 

Figure 2c. Same as in figure 2a for the monochromatic radiation pressure. 

Figure 3. Monochromatic mean intensity at different radii (full lines), together 
with the corresponding blackbody function at T(rj„) (dashed line), for "hot" 
accretion onto a black hole with m = 71. 

Figure 4. Monochromatic mean intensity at different scattering depths (full 
lines) for a "cold", static, plane-parallel atmosphere around a neutron 
star with I = 10~^. The blackbody function at T{Tin) is also drawn for 
comparison (dashed line). 

Figure 5. The emergent spectrum for the model in figure 4 (full line) compared 
with the blackbody at the neutron star effective temperature (dash-dotted 
line) and with the solution obtained by Zampieri et al. (1995, dashed line). 

Figure 6. The gas temperature profile for the model in figure 4 (full line), 
compared with that one found by Zampieri et al. (1995, dashed line). 
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